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ABSTRACT 

We present extended simulations of the large-scale distribution of galaxy clusters in 
several dark matter models, using an optimized version of the truncated Zel'dovich 
approximation (TZA). In order to test the reliability of our simulations, we compare 
them with N-body based cluster simulations. We find that the TZA provides a very 
accurate description of the cluster distribution as long as fluctuations on the cluster 
mass scale are in the mildly non-linear regime (cr 8 < 1). The low computational cost 
of this simulation technique allows us to run a large ensemble of fifty realizations for 
each model, so we are able to quantify accurately the effects of cosmic variance. Six 
different dark matter models are studied in this work: Standard CDM (SCDM), Tilted 
CDM (TCDM) with primordial spectral index n = 0.7, Cold + Hot DM (CHDM) 
with Sl hot = 0.3, low Hubble constant (h = 0.3) CDM (LOWH) and two spatially flat 
low-density CDM models with fl = 0.2 and Aa = 0.8, having two different normaliza- 
tions, cr 8 = 0.8 (ACDMi) and cr 8 = 1.3 (ACDM2). We compare the cluster simulations 
with an extended redshift sample of Abell/ACO clusters, using various statistical mea- 
sures, such as the integral of the two-point correlation function, J3, and the probability 
density function (pdf). We find that the models that best reproduce the clustering of 
the Abell/ACO cluster sample are the CHDM and the ACDMi models. The ACDM 2 
model is too strongly clustered and this is probably overestimated in our simulations 
due to the large cr 8 value of this model. All the other models are ruled out at a high 
confidence level. The pdfs of all models are well approximated by a lognormal distribu- 
tion, consistent with similar findings for Abell/ACO clusters. The low-order moments 
of all the pdfs are found to obey a variance-skewness relation of the form 7 k, S^cr 4 , 
with 53 ~ 1.9, independent of the primordial spectrum shape and consistent with ob- 
servational data. After computing the cluster biasing parameter, b c i, we estimate the 
quantity /3 C ; = fl°' 6 /6 c ; for the different models. Owing to the rather large observational 
uncertainties, (3 C \ = 0.20 ± 0.05, this test does not discriminate strongly between the 
different models. The scale-independence of /3 C ;, and thus of b c i, does, however, sug- 
gest that it is probably a reliable procedure to use the linear biasing model to infer 
the dark matter power-spectrum from observational cluster samples. We also note that 
the abundances of clusters predicted using the Press-Schechter theory provide strong 
constraints on these models: only the CHDM, LOWH and ACDM2 models appear to 
produce the correct number-density of clusters. 



1 INTRODUCTION 

The study of the distribution of matter on the largest scales 
provides important constraints on models of cosmic struc- 
ture formation. If the gravitational instability picture is cor- 
rect, the expected displacements of clusters of galaxies away 
from their primordial positions are much smaller than the 
typical separation of these objects. In principle, therefore, 
clusters of galaxies can yield clues about the primordial spec- 
trum of perturbations that gave rise to them. This is the rea- 
son why so much effort has been devoted to compiling deep 



cluster surveys, starting with the pioneering work of Abell 
(1958), Zwicky et al. (1968) and Abell, Corwin & Olowin 
(1989), and leading up to extended redshift surveys both in 
the optical (e.g. Postman, Huchra & Geller 1992; Dalton et 
al. 1994; Collins et al. 1994, and references therein) and in 
the X-ray (e.g. Nichol, Briel & Henry 1994; Romer et al. 
1994) regions of the spectrum. 

Accompanying the observational challenge of acquiring 
extended cluster redshift surveys, a great deal of effort has 
also been directed towards the provision of reliable statisti- 
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cal characterizations of the cluster distribution. It has been 
established that the cluster two-point correlation function 
is well modelled by a power-law, 

£(,) = (r/r„n. (1) 

Although the slope, 7 ~ 1.8, turns out to be quite similar to 
that of galaxies, the correlation length, r , is much larger (cf. 
Bahcall & Soneira 1983; Klypin & Kopylov 1983; Bahcall 
1988). Different determinations, based on different cluster 
samples, indicate values in the range r = 13-25 h' 1 Mpc* 
(e.g., Nichol et al. 1994, and references therein), while more 
recently the reliability of the power-law model for £(r) has 
also been questioned by different authors (e.g. Olivier et al. 
1993). 

In order to compare the observational cluster data sets 
with different cosmological models several authors have re- 
sorted to large N-body simulations which were designed to 
sample the length-scales relevant to the cluster distribution 
(e.g. White et al. 1987; Bahcall & Cen 1992; Croft & Ef- 
stathiou 1994). The problem with this kind of approach is 
that large N-body simulations are very expensive from a 
computational point of view. Therefore, one is usually forced 
to consider only a limited number of models, with a small 
number of independent realizations for each model. In this 
respect, analytical approaches, based either on Eulerian lin- 
ear theory (e.g. Bardeen et al. 1986; Coles 1989; Lumsden, 
Heavens & Peacock 1989; Borgani 1990; Holtzman & Pri- 
mack 1993) or on the Zel'dovich approximation (ZA, here- 
after; e.g. Doroshkevich & Shandarin 1978; Mann, Heavens 
& Peacock 1993), are in general preferred to numerical sim- 
ulations. Nevertheless, they are of limited utility since one is 
often obliged to resort to oversimplifying assumptions about 
the nature of galaxy clusters. Furthermore, statistics which 
go beyond the two-point correlation function and its Fourier 
transform, the power-spectrum, are hard to handle. Finally, 
it is not clear how shot-noise effects and/or observational 
biases (e.g., redshift-space distortions, selection functions, 
non-trivial sample geometry) can be realistically modelled 
in order to allow a consistent comparison with real data sets. 

In a previous paper (Borgani, Coles & Moscardini 1994, 
hereafter Paper I) we used the truncated Zel'dovich ap- 
proximation (TZA hereafter) to generate cluster simulations 
which were accurate when compared to N-body simulations, 
and, at the same time, computationally so cheap as to en- 
able us to produce many realizations of several dark matter 
(DM) models (see also Blumenthal, Dekel & Primack 1988). 
Recently, Sathyaprakash et al. (1994) compared the TZA 
and several approximations of non-linear gravitational clus- 
tering to direct N-body results. They showed that, although 
the TZA fails to follow small-scale clustering in the multi- 
stream region, it is nevertheless able to account for non-local 
effects due to long wavelength modes in the density fluctu- 
ation spectrum. For this reason, the TZA is not expected 
to provide a correct description of the internal structure 
and mass distribution of non-linear structures like galaxy 
clusters, but it should be accurate in locating their correct 
positions and thus reliably describes their spatial distribu- 
tion. In this respect, it is not necessary to employ the full 



* h is the Hubble constant in units of 100 km s 1 Mpc 1 . 



power and sophistication of modern N-body methods to in- 
vestigate cluster clustering on large scales. It should be clear 
that cluster positions will be accurately determined as long 
as fluctuations on the cluster mass scale are in the linear or 
mildly non-linear regime (os^S 1). However, at later epochs, 
when <78 > 1, shell crossing is no longer negligible on such 
scales and non-linear effects such as cluster infall and merg- 
ing become relevant. 

Given the approximate nature of our method, it is par- 
ticularly important to test its reliability for the task we have 
set it. Indeed, in a recent paper, Gaztafiaga, Croft & Dal- 
ton (1995) have compared cluster simulations based on the 
TZA and on N-body for a model with as = 1.3. Based on 
their analysis of this one model, they reached the general 
conclusion that TZA cluster simulations are not reliable for 
simulating the cluster distribution under any circumstances. 
One of the purposes of the present paper is to demonstrate 
that the objections of Gaztafiaga et al. apply only to mod- 
els in which as is particularly large and that our technique 
is actually extremely accurate for models in which as is of 
order unity or less. To do this, we present cluster simula- 
tions based on an implementation of the ZA that has been 
substantially improved with respect to that used in Paper I 
(see also Plionis et al. 1995, hereafter Paper II). Indeed, we 
can even show that our method is as accurate as iV-body 
methods in situations where its use is appropriate and has 
a fraction of the computational cost. This low cost has al- 
lowed us to generate 50 realizations of each of six different 
initial power-spectra for statistical analysis. This quantita- 
tive analysis is the second purpose of this paper. 

The availability of such a large set of simulations repre- 
sents an extremely powerful test-bed for constraining DM 
models through a detailed comparison of the statistical prop- 
erties of real and simulated cluster distributions. However, 
a reliable quantitative measure of cluster clustering is not 
easy to find. Clusters are rather rare objects with typical 
mean separation of several tens of Mpcs; while bright galax- 
ies have a mean separation comparable to their correlation 
length, clusters have a mean separation which is twice the 
corresponding r value. For this reason, shot-noise effects 
become important on small scales (^ 10 ft -1 Mpc) while, on 
larger scales (^ 40 ft -1 Mpc), a low signal-to-noise ratio is 
expected because the clustering is weak. Robust statistical 
estimators, which are able to provide reliable measures over 
a large range of scales, are required to describe the cluster 
distribution properly and to allow an effective comparison 
with model predictions. 

As a first statistical test, we will use in this paper the 
quantity Js(i?), which is defined through the integral of £(r): 



MR) = 




The advantage of using Js(i?) over £(r) lies in the fact that 
in a sparse distribution of objects, an integral quantity such 
as that defined by eq. (2), should be less susceptible to sta- 
tistical noise than a differential one, such as £(r). 

An alternative method, which is becoming increasingly 
popular, is the study of the probability density function 
(pdf) itself. Usually one attempts to obtain a continuous 
density field by smoothing the discrete distribution of ob- 
jects with some window function (a top-hat or a Gaus- 
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sian are the most commonly used). The smoothing proce- 
dure itself significantly reduces the shot-noise, which could 
dominate the discrete distribution (Gaztafiaga & Yokoyama 
1993). Then one can define the pdf f(g), where q = p/(p), 
and derive its moments, defined by: 



<£"} = J rf(S)dS, (3) 

where 8 = q — \. 

The pdf and moments of different galaxy samples have 
been estimated by various authors (e.g. Saunders et al. 1991; 
Bouchet et al. 1993; Gaztafiaga and Yokoyama 1993; Sheth, 
Mo & Saslaw 1994). Kofman et al. (1994) have compared 
the pdf derived from CDM N-body simulations with that 
of the IRAS sample and the one recovered using the PO- 
TENT procedure with Q = 1 (see also Lahav et al. 1993). 
Their main conclusion is that, if galaxies trace the mass, the 
observed pdf is consistent with Gaussian initial conditions. 

Plionis & Valdarnini (1995, hereafter PV95) studied the 
pdf (and its moments) of the 3-D Abell/ACO smoothed 
cluster distribution and compared them with static simula- 
tions, based on a Gaussian fluctuation spectrum, which re- 
produced the two- and three-point cluster correlation func- 
tions as well as the observed selection effects. They found 
that the real and simulated cluster pdfs are well approxi- 
mated by a lognormal distribution. Cappi & Maurogordato 
(1995) have realized a study of the higher-order moments 
for the discrete Abell/ACO cluster distribution, for both 
projected and redshift samples, while Kolatt, Dekel & Pri- 
mack (1995) estimated the pdf for real cluster samples as 
well as for cluster N-body simulations based both on Gaus- 
sian and non-Gaussian initial CDM fluctuations. They con- 
cluded that no evidence of non-Gaussian initial conditions 
are imprinted into the shape of the cluster pdf. In our Pa- 
per II we compared the variance and the skewness of the 
smoothed Abell/ACO cluster pdf with those obtained from 
the TZA simulations for a list of DM models. In this paper 
we will study the pdf statistics of our latest cluster simu- 
lations. Using the same analysis procedure as that used for 
the Abell/ACO cluster sample (PV95), we will put stringent 
constraints on the models we consider. 

The layout of this paper is as follows. In Section 2 we de- 
scribe the cluster simulations, i.e., how to optimize the ZA, 
the method of cluster identification, the simulation reliabil- 
ity when compared to N-body and higher-order Lagrangian 
results, and the considered models for the power-spectrum. 
In view of the unsuitability of our simulation method for 
studying the small-scale structure of clusters, we also use 
the Press & Schechter (1974; PS hereafter) method to com- 
pare the cluster abundances predicted by DM models with 
the available observational data. This analysis provides an 
independent constraint on the models we are considering. 
In Section 3 we describe the Abell/ACO sample we use. In 
Section 4 we present the analysis of the discrete cluster dis- 
tribution using the J3 integral, while in Section 5 we present 
the pdf and moment analysis of the smoothed cluster dis- 
tribution. In Section 6 we discuss our results and state our 
main conclusions. 



2 THE SIMULATIONS 

2.1 The Zel'dovich approach 

The Zel'dovich approximation (Zel'dovich 1970; Shandarin 
& Zel'dovich 1989) is based on the assumption of laminar 
flow for the motion of a self-gravitating non-relativistic col- 
lisionless fluid. Let q be the initial (Eulerian) position of 
a fluid element and r(q, t) = a(t) x(q, t) the final position 
at the time t, which is related to the comoving Lagrangian 
coordinate x(q, t) through the cosmic expansion factor a{t). 
The ZA amounts to assume the expression 

r(q,i) = a(t) [q + b(t) Vqtf(q)] (4) 

for the Eulerian-to-Lagrangian coordinate mapping. In 
eq.(4) fe(t) is the growing mode for the evolution of linear 
density perturbations and V>(q) is the gravitational poten- 
tial, which is related to the initial density fluctuation field, 
<5(q), through the Poisson equation 

WW - -fg. (.) 

As a result of the factorization of the t- and q- dependence 
in the displacement term of eq.(4), the fluid particles move 
under this approximation along straight lines, with comov- 
ing peculiar velocity 

v(q,i) = x(q,t) = fe(t)VqV>(q). (6) 

Therefore, gravity determines the initial kick to the fluid 
particles through eqs.(5) and (6), and afterwards they do not 
feel any tidal interactions. Particles fall inside gravitational 
wells to form structures, which however quickly evaporate. 
In this sense, the ZA gives a good description of gravita- 
tional dynamics as far as particle trajectories do not inter- 
sect with each other, while its validity breaks down when 
shell-crossing occurs, and local gravity dominates. 

Several prescriptions have been suggested to overcome 
the shortcomings of the ZA, such as adding a small viscous 
term to the equation of motion for the fluid, or by going to 
higher-orders in Lagrangian perturbative theory (e.g. Sahni 
& Coles 1995, and references therein). As a further possi- 
bility, Coles, Melott & Shandarin (1993) have shown that 
filtering out the small-scale wavelength modes in the linear 
power-spectrum reduces the amount of shell-crossing, thus 
improving the performance of the ZA. Melott, Pellman & 
Shandarin (1993) claimed that an optimal filtering proce- 
dure is obtained by convolving the linear power-spectrum 
with the Gaussian filter 

W G (kRf) = e -<^/) 2 / 2 . (7) 

The problem then arises of choosing the filtering radius Rf 
appropriately, in order to suppress shell-crossing as much as 
possible without preventing genuine clustering to build up. 
In Paper I we chose Rf so that the expected mass within 
a Gaussian window of that radius were of the same order 
(~ lO 15 Af0) of the mass for a rich galaxy cluster. The dis- 
advantage of this approach is that it does not rely on any 
objective criterion to optimize the ZA and treats in the same 
fashion different fluctuation spectra, which should produce 
a different amount of shell-crossing. The resulting filtering 
radius, Rf = 5f2 1//3 ft. -1 Mpc, is generally larger than the 
optimal ones, which we use in the present paper, thus caus- 
ing an excessive removal of clustering. 



4 S.Borgani, M.Plionis, P.Coles and L.Moscardini 



Kofman et al. (1994) derived an analytical expression 
for the average number of streams at each Eulerian point, 
N s , as a function of the r.m.s. fluctuation level of the initial 
Gaussian density field. In Figure 1 we plot N s as a function 
of the filtering scale Rf for the six different power-spectra 
that we will consider (see next subsection), evaluated ac- 
cording to eq. (7) of Kofman et al. (1994). As a general 
criterion, we decided to choose Rf for each model so that 
N s = 1.1. We found this to be a reasonable compromise be- 
tween smaller N s values, giving rapidly increasing Rf and 
high suppression of clustering, and larger N s , at which the 
ZA progressively breaks down. The resulting r.m.s. fluctua- 
tion value corresponding to N s = 1.1 is <r = 0.88. 

By adopting this implementation of the TZA, the main 
steps of our cluster simulations are the following: 

(a) Convolve the linear power-spectrum with the Gaussian 
window of eq. (7) and Rf chosen as previously described. 

(b) Generate a random-phase realization of the density field 
on 128 3 grid points for a cubic box of L = 320 ft. -1 Mpc 
aside. 

(c) Move 128 3 particles having initial Lagrangian position 
on the grid, according to the TZA. Each particle carries 
a mass of 4.4 x 10 12 ft _1 Q M . 

(d) Reassign the density and the velocity field on the grid 
through a TSC interpolation scheme (e.g. Hockney & 
Eastwood 1981) for the mass and the moment carried by 
each particle. 

(e) Select clusters as local density maxima on the grid ac- 
cording to the following prescription. If d c i is the average 
cluster separation, then we select N c i = (L/d c i) 3 clusters 
as the N c i highest density peaks. In the following, we as- 
sume d c i = 40 h -1 Mpc, which is appropriate for the com- 
bined Abell/ACO cluster sample to which we will compare 
our simulation results (see Section 3). Therefore, we will 
analyze a distribution of 512 clusters in each simulation 
box, with periodic boundary conditions. 



2.2 Reliability of the TZA 

Before entering into the presentation of our analysis, we 
present the results of detailed tests of the reliability of TZA 
for simulating the large-scale distribution of galaxy clusters. 
For this test, we ran a PM simulation for an initial spectrum 
corresponding to the Cold+Hot DM (CHDM) model with 
Qhot = 0.3 (see Section 2.3 for more details). The size of 
the simulation box and the number of grid points and parti- 
cles were identical to those for the TZA simulations. We do 
not attempt here to distinguish between hot and cold parti- 
cles, since the simulation is intended only for a comparison 
with the Zel'dovich approach. In any case, it is reasonable 
to assume that the adopted mass resolution is low enough 
for effects of neutrino free-streaming to be negligible. Clus- 
ters are identified as local maxima on the grid, following the 
same prescription outlined in Section 2.1. We compare the 
outputs of the N-body simulation to those of several TZA 
realizations, each based on the same initial phase assign- 
ment, but having different filtering radii Rf. 

In Figure 2 we compare the cluster distributions within 
a slice 80 ft. -1 Mpc thick and the cluster count -in-cell vari- 
ance, <j 2 (i?), for PM and TZA simulations, at two differ- 
ent evolutionary stages, corresponding to as = 0.67 and 



<J8 = 1 for the r.m.s. fluctuation amplitude within a top-hat 
sphere of 8 ft. -1 Mpc radius. The less-evolved stage is con- 
sistent with the two-year COBE normalization supplied by 
Bennett et al. (1994), but the later stage corresponds to a 
higher normalization than this. We consider the second stage 
only in order to assess the reliability of the TZA simulations 
when non-linear effects appear on the cluster mass scale. 
The variance has been estimated according to 
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where (-^V}^ and (-^V 2 }^ are the average and the second 
order moment, respectively, for counts within 20,000 ran- 
domly placed spheres of radius R and the second term in 
the l.h.s. of eq.(8) represents the correction for Poissonian 
shot-noise. We take R between ~ 10 and 80 ft." 1 Mpc, since 
smaller scales are heavily affected by shot-noise and larger 
scales by effects of periodic boundary conditions. 

At as = 0.67 the TZA reproduces the variance of the 
PM cluster distribution remarkably well, especially when 
the linear power-spectrum is mildly filtered with Rf = 
3 ft -1 Mpc, a value which is quite close to that required for 
N s = 1.1. For larger filtering no allowance is made for gen- 
uine clustering to build up. As a consequence, the large- 
to-small scale power ratio is decreased, with a subsequent 
increase of the cluster correlation. On the other hand, a sim- 
ilar clustering suppression is caused by shell-crossing when 
a smaller filtering, or no filtering at all, is taken. It should be 
stressed, however, that the agreement between N-body and 
TZA cluster simulations is not merely in a statistical sense. 
In the upper panel we superimpose the cluster distributions 
from PM (filled dots) and TZA (R f = 3 ft _1 Mpc; open cir- 
cles) simulations. The relative positions are identical for al- 
most all the clusters. This confirms that the TZA places 
peaks in the correct positions when it is suitably "trun- 
cated", thus disproving the general conclusion by Gaztariaga 
et al. (1995) that the TZA is not reliable enough for cluster 
simulations. 

In order to verify that the close agreement between 
TZA and PM cluster simulations is generally attained, even 
when considering N-body simulations with different dynam- 
ical and mass resolution, we compare in Figure 3 the two- 
point cluster correlation function for our CHDM simulations 
(open dots) to that obtained by Klypin & Rhee 1994; KR94, 
hereafter) by evolving until as = 0.67 the same spectrum 
with a PM N-body code with 256 3 grid points and particles 
within a box of 200 ft" 1 Mpc. The KR94 results are obtained 
as an average over 2 realizations and error bars are quasi- 
Poissonian sampling uncertainties while our results refer to 
the average taken over 50 random realizations. Our error 
bars are estimated as the r.m.s. scatter over this ensemble. 
By comparing this plot with Fig. 2 of Paper I, it is appar- 
ent that we have improved our implementation of the TZA 
by increasing the resolution and by optimizing the power- 
spectrum filtering. The agreement between the TZA and 
KR94 results is really remarkable and further confirms the 
reliability of the TZA to follow correctly the mildly non- 
linear clustering regime. 

This agreement, however, is not so good when later evo- 
lutionary stages are considered. Indeed, at as = 1, although 
<j 2 (i?) for TZA clusters increases to some extent, it is rather 
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stable or marginally decreasing for PM clusters. As a re- 
sult, the TZA systematically overestimates the cluster cor- 
relations for any choice of the filtering radius. The lack of 
evolution of a 2 (R) in the N-body simulations is more likely 
to be ascribed to a balance between two competitive effects, 
the infall of close clusters, which tends to enhance the small- 
scale clustering, and cluster merging, which decreases the 
number of close pairs, rather than to a freezing of the clus- 
tering. This kind of picture is also supported by studies of 
the correlation function of cluster peculiar velocities (Cen, 
Bahcall & Gramann 1994; Croft & Efstathiou 1995), which 
has negative values below 20 ft. -1 Mpc, showing therefore ev- 
idence of infall on such scales. In the N-body treatment and 
in terms of the cluster distribution, the disappearance of 
close pairs within overdense regions leads one to pick up rel- 
atively more clusters in the "field", as one is willing to select 
a fixed number of such clusters (cf. the upper right panel in 
Figure 2). 

We conclude from this analysis that the TZA is very 
accurate for cluster simulations as long as effects of non- 
linear gravitational dynamics are negligible on the cluster 
mass scale (i.e., os^S 1). Although this holds in most cases 
of cosmological relevance (for instance, as — 0.6 is required 
for n o = 1 CDM like models to fit the observed cluster 
abundance; cf. White, Efstathiou & Frenk 1993, and Sec- 
tion 2.4 below), there are nevertheless a few exceptions, for 
which the conditions for our method to be accurate are not 
satisfied: we should therefore bear in mind the above words 
of caution. 

On the basis of this result, one may wonder whether in- 
cluding gravitational effects in the ZA, at least to first order, 
would improve the agreement with N-body for as = 1. To 
this purpose, we realised simulations based on the second- 
order Lagrangian theory (e.g., Moutarde et al. 1991; Bouchet 
et al. 1992; Buchert & Ehlers 1993; Gramann 1993; Cate- 
lan 1995). Expanding equation (4), for the particle displace- 
ments, to the second order in fe(t) we obtain 

r (q, t) = a(t) [q + 6(i)p(q) + b 2 (t)*(q)] . (9) 

Here, p(q) = VqV>(q) is the first order term of eq.(4), while, 
for n = 1, the term 

l<t<j<3 

includes tidal effects, which generate the first-order devia- 
tion from the straight-line particle trajectories. 

In Figure 4 we plot the cluster two-point correlation 
function, £(r), for the CHDM simulations, based on the 
first-order (filled circles) and on the second-order (open 
circles) TZA, for both as = 0.67 (left panel) and as = 1 
(right panel). Results correspond to the average over 10 re- 
alizations and error bars are la scatter over this ensemble. 
Note that going to second-order in Lagrangian theory does 
not introduce any significant change in the cluster distribu- 
tion. This suggests that non-linear gravitational effects, like 
merging and infall of structures, cannot be treated with a 
perturbative approach, at least at this order. 



2.3 Dark matter models 

We ran simulations for six different models of the initial fluc- 
tuation spectrum. For each model, we generate 50 random 
realizations, so as to reliably estimate the effect of cosmic 
variance. All the models, except the ACDMi one (see be- 
low), are normalized to be consistent with the COBE mea- 
sured quadrupole of CMB temperature anisotropy (Bennett 
et al. 1994). The models we have considered are the follow- 
ing. 

(1) The standard CDM model (SCDM), with U = l,h = 0.5 
and as = 1 for the r.m.s. fluctuation amplitude within a 
top-hat sphere of 8 h -1 Mpc. The rather large normaliza- 
tion of this model could produce an overestimate of cluster 
clustering in the light of the previous discussion about the 
reliability of the TZA cluster simulations. However, as we 
shall see, this model already produces too weak clustering 
for clusters, while a more accurate treatment would even- 
tually produce even weaker correlations. Therefore we will 
give at most only an underestimate of the confidence level 
to which SCDM should be ruled out. 

(2) A tilted CDM model (TCDM), with n = 0.7 for the 
primordial spectral index. Tilting the primordial spectral 
shape from the scale-free one has been suggested in order 
to improve the CDM description of the large-scale struc- 
ture (e.g. Cen et al. 1992; Tormen et al. 1993; Liddle & 
Lyth 1993; Adams et al. 1993; Moscardini et al. 1995). 

(3) A low Hubble constant CDM model (LOWH), with h = 
0.3. Decreasing the Hubble constant has the effect of in- 
creasing the horizon size at the equivalence epoch, thus 
pushing the turnover of the spectrum to its scale-free form 
out to larger scales. The relevance of this model in alle- 
viating several cosmological problems has been recently 
emphasized by Bartlett et al. (1994). 

(4) A Cold + Hot DM model (CHDM), with U hot = 0.3 
for the fractional density contributed by the hot particles. 
For a fixed large-scale normalization, adding a hot com- 
ponent has the effect of suppressing the power-spectrum 
amplitude at small wavelengths (see, e.g., Klypin et al. 
1993; Klypin, Nolthenius & Primack 1995, and references 
therein, for the relevance of CHDM). Although the small- 
scale peculiar velocities are lowered to an adequate level, 
the corresponding galaxy formation time is delayed so 
that such a model is strongly constrained by the detec- 
tion of high-redshift objects (e.g. Ma & Bertschinger 1994; 
Klypin et al. 1995, and references therein). 

(5) A spatially flat, low-density CDM model (ACDMi), with 
Q = 0.2, Qa = 0.8 for the cosmological constant term 
(e.g., Bahcall & Cen 1992; Baugh & Efstathiou 1993; Pea- 
cock & Dodds 1994) and as = 0.8. With this normal- 
ization, ACDMi has a significantly lower amplitude than 
that implied by COBE data. 

(6) The same model as in (5), but with a larger normal- 
ization, as = 1.3 (ACDM2), so as to be consistent 
with COBE results. Having in mind the comparison be- 
tween PM and TZA simulations previously discussed, the 
ACDM2 model is expected to have a too large as value to 
be adequately treated by the TZA. For this reason, we pre- 
fer to consider also the ACDMi model, whose lower nor- 
malization allows it to be properly handled by the TZA. 
Since cluster correlations have been shown to be almost 
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Table 1. The models. Column 2: the density parameter Qq ; Col- 
umn 3: the cosmological constant term ; Column 4: the density 
parameter of the hot component Clhot\ Column 5: the primordial 
spectral index n; Column 6: the Hubble parameter h; Column 7: 
the linear r.m.s. fluctuation amplitude at 8 h~ 1 Mpc as ; Column 
8: The filtering radius, Rf, in units of /i — 1 Mpc, corresponding to 
N s = 1.1 for the level of orbit crossing. 



Model 


f2o 






n 


h 


<T8 


Rf 


SCDM 


1.0 


0.0 


0.0 


1.0 


0.5 


1.0 


4.4 


TCDM 


1.0 


0.0 


0.0 


0.7 


0.5 


0.5 


1.6 


LOWH 


1.0 


0.0 


0.0 


1.0 


0.3 


0.6 


2.4 


CHDM 


1.0 


0.0 


0.3 


1.0 


0.5 


0.7 


2.2 


ACDMi 


0.2 


0.8 


0.0 


1.0 


1.0 


0.8 


3.2 


ACDM 2 


0.2 


0.8 


0.0 


1.0 


1.0 


1.3 


6.3 



independent of as, we could argue that results at as = 0.8 
can be considered as representative of those at as = 1.3. 

The transfer functions for the above models have been 
taken from Holtzman (1989), except that of LOWH, which 
is taken from Bond & Efstathiou (1984), with suitably cho- 
sen shape parameter T = Q h = 0.3. We note that the 
latter transfer function assumes the baryonic component to 
be negligible. However, for h = 0.3, constraints from the Big 
Bang nucleosynthesis imply a considerable baryon fraction, 
U b ~ 0.15 (e.g., Walker et al. 1991). Although this could be 
a limitation for small scale considerations, it hardly affect 
the results at the large scales we are interested in. All the 
model parameters are listed in Table 1. 

In Paper II we also considered a purely open CDM 
model with Q = 0.2, without cosmological constant and 
with as = 1. However the only difference between low- 
density models with and without cosmological constant lies 
in the different way one normalizes to the large-scales CMB 
anisotropics, while the spectrum shapes are virtually identi- 
cal. Therefore, we expect this model to be just intermediate 
between the two above ACDM models. 

Each power-spectrum is suitably smoothed on the scale 
Rf according to the prescription described in Section 2.1. 
In Figure 1 we plot the average stream number per Eulerian 
point, N s , as a function of the filtering radius for all the 
above models. The intersection of the N s = 1.1 line with 
each curve indicates the smoothing scale adopted for the 
corresponding spectrum (see also Table 1). Note that the 
larger as value for ACDM2 requires a stronger filtering to 
suppress shell-crossing. 

In Figure 5 we plot the projected particle distribution 
within a slice 10 ft. -1 Mpc thick, superimposing the cluster 
distribution, for the CHDM model. In order to better show 
how the identified clusters trace the underlying density field, 
we used simulations within a 640 h -1 Mpc box. As expected, 
clusters are preferentially located at the knots corresponding 
to the intersections between filaments, while they avoid long 
filaments and flattened pancakes, so as to generate devoid 
regions of size ~ 200 h -1 Mpc, at least within the slice. 



2.4 Cluster Abundances 

As an independent constraint on the above models, we have 
computed the expected cluster abundances, as predicted 
by the standard PS formalism. If structures are identified 
through a filter W on a scale R, in order to have a mass 
M = fpR 3 (p is the average matter density), the PS for- 
mula for the number density of objects with mass between 
M and M + dM is 



n(M)dM 
where 

v{R) = 

a\R) = 



1* C< 

* Jr 



v(R) 

a(R) 



exp 



2<t 2 (R) 



(11) 
R 2 ' ( ' 



2^a\R) j 
2^f 



k* P(k) 



dW 2 {kR) dk_ 
d(kR) kR 



k 2 P(k)W 2 (kR)dk. 



(12) 



Therefore, the total abundance of objects of mass larger than 
M is 



N(> M) 



p CO 

/ n(M')dM'. 

J M 



(13) 



In the above expressions / is a "form factor", which depends 
on the shape of the filter W : f = (27r) 3/ ' 2 for a Gaussian fil- 
ter, and / = 47r/3 for a top-hat filter. A Gaussian filter will 
be assumed in the following analysis. The parameter S c is 
the critical density contrast, which represents the threshold 
value for a fluctuation to turn into an observable object, 
if evolved to the present time by linear theory. Arguments 
based on a simple spherical collapse suggest S c = 1.68, but 
the inclusion of non-linear effects, as well as aspherical col- 
lapse, may lead to a lower value of S c . For example, Klypin 
& Rhee (1994; KR94 hereafter) found that the cluster mass 
function in their CHDM N-body simulations is well fitted 
by eq.(ll) by taking S c = 1.5. 

White, Efstathiou & Frenk (1993) resorted to X- 
ray data for the temperatures of the gas component 
of clusters and estimated a cluster abundance of about 
4 x 10 _6 ( ft. _1 Mpc) -3 for masses exceeding M = 4.2 x 
10 h~ Mq. Using observed cluster velocity dispersion, Bi- 
viano et al. (1993) obtained an abundance of about 6 x 
10 -6 ( ft. -1 Mpc)~ 3 for clusters exceeding the above mass 
limit . 

In Figure 6 we compare model predictions at different 
S c values to the above observational estimates. Note that 
realistic uncertainties on cluster abundances are probably 
larger than the difference between the two above values. 
They should include variations in the average cluster num- 
ber density between different samples, biases toward high 
mass for observations of cluster velocity dispersion, uncer- 
tainties in the model used to relate gas temperature and 
cluster mass, etc. Keeping in mind such warnings, we note 
from Figure 6 that ACDMi is ruled out; for S c = 1.5 it 
produces more than one order of magnitude less clusters 
than observed. This agrees with the suggestion of White 
et al. (1993), that a higher normalization [as — 1.4) is re- 
quired for such a low-density model to provide a correct 
number of clusters. Consistently, ACDM2 gives the right 
cluster abundance for reasonable values of S c . As for the 
Q = 1 models, it turns out that the resulting abundances 
depend mostly on the as normalization value and not on 
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the shape of the spectrum. This is not surprising: the Gaus- 
sian smoothing scale, R ~ 4.6 ft. -1 Mpc, which encompasses 
a mass of 4.2 x 10 14 ft. -1 Mq, is equivalent to a top-hat sphere 
of about 7 h -1 Mpc, which is rather close to the normaliza- 
tion scale. As a result, SCDM turns out to produce too many 
clusters for any reasonable value of S c . On the other hand, 
the low normalization of TCDM turns into a severe under- 
production of clusters, even at the smallest S c values. The 
only two flat models which generate cluster abundances in 
agreement with the results by White et al. (1993) and Bi- 
viano et al. (1993) for a reasonable choice of S c are LOWH 
and CHDM. 

A note of caution about the strength of the constraints 
emerging from the Press-Schechter analysis, is due for a 
number of reasons. First, there is some evidence of a dis- 
crepancy between the mass profiles of clusters inferred from 
X-ray data (e.g. Edge & Stewart 1991) and from gravita- 
tional lensing considerations (e.g. Kaiser et al. 1994), so it 
is not clear whether the cluster mass function inferred from 
the X-ray data is correct. Likewise, there is a possibility 
that the distribution of cluster peculiar velocities may be 
affected by subclustering. One should also mention that the 
applicability of the Press-Schechter method is itself open to 
some doubt. Although, as we mentioned above, it appears 
to perform well for the CHDM model when compared with 
N-body simulations, its accuracy is yet to be verified for the 
other models. In the case of LOWH, where a much higher 
fraction of the cluster mass is baryonic than in the other 
models, one might imagine this formalism to be particularly 
suspect. We therefore take the constraints emerging from 
this analysis to be indicative but not watertight. 



3 THE CLUSTER SAMPLE 

We use the combined Abell/ACO R > cluster sample, as 
defined in Plionis & Valdarnini (1991) [hereafter PV91] and 
analysed in Plionis, Valdarnini & Jing (1992) [hereafter PVJ] 
and Plionis & Valdarnini (1995). The northern sample, with 
dec > —17° (Abell), is defined by those clusters that have 
measured redshift z^S 0.1, while the southern sample (ACO; 
Abell, Corwin & Olowin 1989), with dec< -17°, is defined 
by those clusters with mio < 16.4 (note that with this defini- 
tion and due to the availability of many new cluster redshifts 
only 7 ACO clusters have mio estimated redshifts from the 
mio — z relation derived in PV91). Both samples are limited 
in Galactic latitude by \b\ > 30°. The redshifts have been 
taken from a number of studies, the references of which can 
be found in PV95. The total number of clusters in our sam- 
ples is 357 and 157, for Abell and ACO respectively. 

To take into account the effect of Galactic absorption, 
we assume the usual cosecant law: 

P(\b\) = dex [a(l - csc|6|)] (14) 

with a PS 0.3 for the Abell sample (Bahcall & Soneira 1983; 
Postman et al. 1989) and a PS 0.2 for the ACO sample (Ba- 
tuski et al. 1989). The cluster-redshift selection function, 
P(z), is determined in the usual way (cf. Postman et al. 
1989; PVJ; PV95), by fitting the cluster density, as a func- 
tion of z (see the above reference for details). Cluster dis- 
tances are estimated using the standard relation: 



R = H oq il + z) + V ~ "X 1 " v^TT)] («) 

with H = 100 ft km sec" 1 Mpc" 1 and q = o /2. Strictly 
speaking, eq.(15) holds only for vanishing cosmological con- 
stant. Therefore, for a consistent comparison with the sim- 
ulation models, we should use different R-z relations for 
the Abell/ACO analysis. However, we verified that final re- 
sults are essentially independent of the choice of the (A, Q ) 
parameters used in the simulations. For this reason, in the 
following we will present results for real data only based on 
assuming eq.(15) with q = 0.2. 

PVJ and PV95 found that the Abell and ACO clus- 
ter number densities, out to their limit of completeness, 
are ~ 1.4 x 10" 5 ft 3 Mpc" 3 and ~ 2.1 x 10" 5 ft 3 Mpc" 3 , 
corresponding to mean separations d c i ss 41 ft" 1 Mpc and 
d c i PS 36 ft" 1 Mpc, respectively. The higher space-density of 
ACO clusters is partly due to the huge Shapley concentra- 
tion (Shapley 1930), but a significant part is also due to 
systematic density differences between the Abell and ACO 
cluster samples, as a function of z, which has been noted 
in a number of studies (cf. PV91 and references therein) 
and which could be attributed to the high sensitivity of the 
Illa-J emulsion plates. In PV95, this effect was taken into 
account by normalizing the densities of the two samples us- 
ing a radial matching function, W(R), which is defined as 
the ratio between the average densities for Abell and ACO 
clusters at equal volume shells. 

In the following, we compare results based on the 
Abell/ACO sample with those derived from our simulated 
cluster populations, selected so that d c i = 40 ft" 1 Mpc. Vari- 
ations in d c i of the order of the Abell-ACO difference, does 
not significantly affect the resulting statistical properties. Fi- 
nally, following PV95 we restrict our analysis of the real clus- 
ter sample within a maximum distance of R ma x = 240 ft" 1 
Mpc, in order to minimize the uncertainties due to the ap- 
proximate character of the redshift selection function, P(z), 
and of the radial matching function, W(R), especially at 
large distances. 



4 STATISTICS OF THE DISCRETE CLUSTER 
DISTRIBUTION 

Our first statistical test for comparing the real data and the 
simulations involves the evaluation of the quantity Js(i?), 
defined by eq. (2). It is straightforward from the definition 
of this quantity to construct the estimator 

^> = t(t-')' (16) 

where N n b is the average number of cluster neighbours 
within a distance R from a cluster, while N is the expected 
number of neighbours for a random cluster distribution (es- 
timated at the positions of the real clusters). Therefore, 
Jz(R) cx R 3 ' as long as £(r) cx r '. 

It has been argued (cf. KR94) that the scale at which 
the power-law shape of £(r) breaks and firstly crosses zero, 
is a potentially powerful test for cosmological models. How- 
ever, since such a scale corresponds by definition to the weak 
clustering regime, its detection can be heavily affected by 
statistical noise; for an explicit demonstration of this, see 
Paper I. In this respect, the analysis of J3 should have the 
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Table 2. J^^R) values at different scales (in h 1 Mpc units) for 
both simulated and real cluster distributions. 



Model 


M R ) 


x lQ-^h^Mpc) 3 






28.9 


43.2 


64.3 




SCDM 


5.00 ± 0.64 


5.55 ± 1.34 


3.86 ± 2 


12 


TCDM 


4.93 ± 0.74 


6.05 ± 1.60 


5.71 ± 2 


74 


LOWH 


5.12 ± 0.79 


6.10 ± 1.60 


5.16 ± 2 


83 


CHDM 


8.39 ± 0.94 


10.19 ± 1.98 


8.77 ± 3 


36 


ACDMi 


8.47 ± 1.06 


11.04 ± 2.23 


11.67 ± 4 


33 


ACDM 2 


18.41 ± 1.96 


23.88 ± 3.03 


26.36 ± 6 


91 


Abell/ACO 


8.10 ± 1.69 


10.28 ± 2.53 


6.96 ± 5 


17 



to be an overestimate, in the light of the considerations dis- 
cussed in Section 2.2. Therefore, if we follow the suggestion 
that cluster correlations are nearly independent of as (cf. 
Croft & Efstathiou 1994; see also Section 2.2), and we take 
the results of ACDMi at as = 0.8 to be representative of 
those at as = 1.3, we can conclude that the ACDM model 
is viable as far as cluster clustering is concerned (see also 
Dalton et al. 1994; Gaztanaga et al. 1995). Note that the 
results for TCDM and LOWH are remarkably similar. This 
agrees with the expectation that, as far as the shape of the 
power-spectrum is concerned, a change in the Hubble pa- 
rameter ft is roughly equivalent to a change in the spectral 
index n according to the relation Aft. = —An (cf. Lyth & 
Liddle 1994). 



advantage of being more stable and suffering less from ob- 
servational biases. 

In Figure 7 we plot Js(i?) for the real data (filled cir- 
cles) and simulations (open circles). We estimated N for 
the real data by averaging over 100 random samples, having 
the same selection criteria (boundaries, galactic extinction 
function, redshift selection and systematic Abell/ACO dif- 
ferences) as the real one. Error bars for the simulated sam- 
ples are la scatter over the ensemble of 50 realizations. In 
Table 2 we report values of J3 for data and simulations at 
three different scales. The quoted uncertainties for real clus- 
ter analysis are la scatter estimated over an ensemble of 
100 bootstrap resamplings. Such errors are not plotted in 
Figure 7. In fact, since we are asking which is the proba- 
bility that a given model generates a result like that of the 
observed cluster distribution, its 'success' is just measured 
by the distance of the real data point from the cosmic r.m.s. 
error bars. This should be taken into account when judg- 
ing to which confidence level a model has to be accepted or 
rejected. In Paper II we verified that intermediate scales of 
few tens of Mpcs are best suited to constrain DM models, 
when using the cluster distribution, smaller and larger scales 
being affected by shot-noise and low signal-to-noise ratio, 
respectively. Furthermore, richness contamination of cluster 
correlations should have effects only on rather small scales, 
below 10 h~ 1 Mpc (see Olivier et al. 1993). For these rea- 
sons, in the present analysis we do not consider scales much 
smaller than 20 ft -1 Mpc as well as larger than 60 ft -1 Mpc. 

For the Abell/ACO sample, Js(i?) increases up to 
R ~ 35 ft -1 Mpc, flattens at a scale corresponding to the 
break of the power-law shape of £(r), and eventually de- 
clines at i?^ 50 ft -1 Mpc, after which £(r) becomes nega- 
tive. By comparing this result with those of the simulations, 
it turns out that the only models which overcome this test 
are CHDM and ACDMi , although both of them seem to 
produce too strong clustering at smal scale and the second 
gives a marginal clustering excess at the largest scale. All the 
other models are ruled out at a > 2a level. The cluster dis- 
tributions for the SCDM, LOWH and TCDM models are too 
weakly clustered over the whole scale range. Note that the 
SCDM has a rather flat Js(i?) profile, according to the ex- 
pectation that this model has a cluster two-point correlation 
function which declines rapidly beyond ~ 20 ft -1 Mpc. Con- 
versely, ACDM2 generates too much clustering, with Js(i?) 
increasing up to i?^ 60 ft -1 Mpc. However, we should keep in 
mind that the large clustering detected in this model ought 



5 STATISTICS OF THE SMOOTHED 
CLUSTER DISTRIBUTION 

In order to facilitate the comparison of our results with those 
obtained by PV95 from the combined Abell/ACO cluster 
sample, we followed basically the same procedure as they 
did, and which we briefly describe below. We obtain a con- 
tinuous cluster density field by smoothing the cluster distri- 
bution on a grid, with grid-cell width of 20 ft -1 Mpc (16 3 
grid-points), using a Gaussian kernel: 



W(|x i -x 9 |)= (2*R 2 m ) 3/2 exp 



x; 



(17) 



The smoothed cluster density, at the grid-cell positions x 9 , 
is then: 



p( x 9 ) 



£. p(x,)>V(|x,-x g |) 

/Ml' 



ig\)d 3 X 



(18) 



where the sum is over the distribution of clusters with po- 
sitions x;. In order to study the cluster density field at dif- 
ferent smoothing scales, we use three radii for the Gaussian 
kernel: R sm = 20, 30 and 40 ft" 1 Mpc with |x;-x 9 | < 3i? sm . 
Therefore the integral in the denominator of eq. (18) has a 
value smaller than unity (~ 0.97). 



5.1 The probability density function 

As a first test for the smoothed cluster density field, we 
work out the probability density function, f(g), which rep- 
resents a low-order (one-point) statistics. We then compare 
the pdf of each set of cluster simulations with the observed 
Abell/ACO pdf, derived by PV95, as well as with the fol- 
lowing theoretical models. 



(a) The Gaussian distribution given by 



fie) 



exp 



2a 2 



(19) 



where a is the standard deviation of q (= p/(p)). If f(g) is 
a Gaussian then it should be defined in an infinite interval, 
which implies that f{g < 0) ^ 0. Since, however q > 
by definition, f(g) is expected to be well approximated 
by a Gaussian only in the limit a — > 0. In this case the 
skewness, 7 (= (<5 3 )), vanishes. Even in the case of an 
initial Gaussian density field, the gravitational evolution 
acts in such a way as to increase the variance a 2 , and 
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thus, due to the constraint q > 0, f(g) has to become 
positively skewed. For as long as the variance a 2 = {S 2 } is 
small, the deviation of the pdf shape from a Gaussian is 
well approximated by the Edgeworth expansion (Colombi 
1994). 

(b) The lognormal distribution given by 



fie) 



1 



: exp 



(lng - n L ) 2 



2<t| 



(20) 



where q is obtained through an exponential transforma- 
tion of a Gaussian random variable % as q = exp(%). In 
eq. (20), fiL and ox are the mean and standard deviation 
of In q respectively. It has been argued that this distri- 
bution describes the distribution of density perturbations 
resulting from Gaussian initial conditions in the weakly 
non-linear regime Coles & Jones 1991). Bernardeau & 
Kofman (1994) have shown that the lognormal distribu- 
tion is not really a natural consequence of mildly non- 
linear gravitational evolution, but a very convenient fit 
only in some portion of the (<r, ra)-plane (i.e. <r <C 1 and 
spectral index n PS — 1). It has nevertheless been found to 
give an extremely good fit to the CDM density and the 
IRAS galaxy pdf in the weakly-linear regime (Kofman et 
al. 1994), as well as to the observed Abell/ACO cluster 
distribution (PV95). 

The pdf resulting from the application of the Z A to Gaus- 
sian initial fluctuations (Kofman et al. 1994): 

9 x 5 3 ' 2 



fie) 



p CO 



ds e 



-(s-3) 2 /2<r : 



-Pi /2a 



+ e~ 



\/2a 



\/2o 



/3 n (s) = -\/5 s {1/2 + cos[2/3 (n - 1)tt 
+ 1/3 arccos (54/gs 3 - l)]} , 



(21) 



where a is the rms amplitude of density fluctuations and 
N s is the average stream number per Eulerian point. 

Note that discreteness effects could be important since 
they affect the shape of the pdf and the estimation of its mo- 
ments; especially at small R sm , when the number of clusters 
in the Gaussian sphere is small, and/or when the smoothing 
fails to create a continuous density field due to discreteness 
(in our case this is apparent in the i? sm = 20 ft. -1 Mpc case 
for q < 0.8). In the case of a Poisson sampling of an underly- 
ing continuous density field, the shot-noise contributions to 
the moments can be easily estimated and corrected for (cf. 
Peebles 1980). The cluster distribution cannot, however, be 
meaningfully regarded as a Poisson sampling of the under- 
lying (galaxy) distribution, since clusters are expected to 
form only at high density peaks. Consequently the Poisso- 
nian shot-noise correction could not give a reasonable de- 
scription of discreteness effects (Coles & Frenk 1991; Bor- 
gani et al. 1994). Moreover, Gaztafiaga & Yokoyama (1993) 
have shown that the smoothing process itself considerably 
suppresses these shot-noise effects. For these reasons PV95 
did not use any shot-noise corrections. To make a consistent 
comparison of our models with the data, we also did not in- 
clude such corrections in our analysis. Since all the model 
cluster distributions have the same mean number density 



and we treat them similarly, the possible effects of shot- 
noise are accounted for in the same way in both the data 
and the simulations: we are therefore comparing like with 
like. 

In Table 3 we present the results of the comparison be- 
tween the simulation pdf, the PV95 Abell/ACO cluster pdf, 
the lognormal distributions and a Gaussian distribution, us- 
ing a x 2 



E 



/sim( \ rtheor ( \ 

i je) - fi je) 



(22) 



where the weights e 2 correspond to cosmic variance. Note 
that we derived the simulation pdf in redshift space so that 
a consistent comparison with the PV95 results can be made; 
the comparison with theoretical models is done in real space. 

In Figure 8 we present the simulation cluster pdfs for 
the CHDM and SCDM models together with the PV95 
Abell/ACO pdf at R sm = 20 and 40 ft." 1 Mpc, and in Figure 
9 we make the comparison between these two cosmologi- 
cal models and the theoretical distributions. The error bars 
represent the scatter around the ensemble mean values (cos- 
mic variance). As in Table 3, comparisons with theoretical 
models and real data are made in real space and in redshift 
space, respectively. There is an excellent agreement between 
the CHDM and the Abell/ACO cluster pdfs while in the 
SCDM case there is a clear discrepancy at small (< 0.5) 
and large £>'s. 

It is apparent that: 

(i) The Gaussian distribution of does not provide a good 
fit at any i? sm < 30 ft -1 Mpc and for any model (cf. 
column 4 of Table 3). For R sm = 40 ft" 1 Mpc the Gaus- 
sian fit is acceptable only for the SCDM, LOWH and 
ACDM 2 models. 

(ii) The Zel'dovich pdf model of eq.(21) is inconsistent with 
the simulation results, even though the underlying dy- 
namics governing the cluster distribution are described 
by the ZA. We find that this distribution is ruled out 
at a confidence level larger than 99.99%, for R sm < 30 
ft -1 Mpc for all the simulation models. This is the rea- 
son why we did not show results for this model in Table 
3. The SCDM, TCDM, CHDM and ACDMi models are 
only consistent at a ~ 20% — 25% level for R sm = 40 
ft -1 Mpc. One may argue that, since the ZA pdf is de- 
signed to describe the DM clustering, its failure for the 
cluster distribution is nothing but the consequence of 
not accounting for the mass within clusters in the anal- 
ysis of their distribution. To check this, we repeated 
the analysis by weighting simulated clusters according 
to their mass and found no appreciable differences in 
the pdf shapes. 

(iii) All simulation pdfs are well approximated by lognor- 
mal distribution of eq.(20), irrespective of their differ- 
ent power-spectra, especially when the larger smooth- 
ing radii are considered (cf. column 3 of Table 3). Note 
that at the i? sm = 20 ft -1 Mpc case the comparison 
is done for q > 0.8 because, at lower values of q, dis- 
creteness effects introduce significant noise. Therefore, 
in contrast to the case of the matter distribution, the 
lognormal fit to the cluster pdf is more likely to be con- 
nected with the high-peak biasing description of cluster 
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Table 3. % 2 probabilities that the indicated simulation model 
pdf could have been drawn from a parent distribution given by 
the lognormal, Gaussian or the real cluster (PV95) pdf. For the 
R srn = 20 h~ 1 Mpc case the comparison is performed for q > 0.8 
(see text). Also reported in Column 6 are the values of the reduced 
skewness S3 for simulations and for real data. 



Model 


R sm 


P 


LN 
K 2 


pG 

r x 2 








5 3 








20 





02 





00 





01 


1 


87 


± 





29 


SCDM 


30 





83 





00 





75 


1 


83 


± 





70 




40 





99 





92 





97 


1 


60 


± 


1 


25 




20 





15 





00 





00 


1 


96 


± 





32 


TCDM 


30 





99 





00 
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formation than being due to non-linear gravitational ef- 
fects, which dominate much smaller scales (Bernardeau 
& Kofman 1994). 

The scale which best discriminates between different 
models and the Abell/ACO data is clearly R sm = 
20 ft. -1 Mpc (cf. column 5 of Table 3). At larger scales, 
all the models, except ACDM2, produce acceptable fits. 
The only models that produce a pdf consistent with the 
PV95 results, at all 3 smoothing radii, are CHDM and 
ACDMi , with the former performing however system- 
atically better. 



5.2 Moments of the pdf 

According to eq. (3), the moments of the pdf give a large 
weight to the high density tail (<5 > 1) of the pdf. They 
are therefore expected to suffer less from shot-noise effects, 
which are smaller in the overdense parts of the distribution. 
In Paper II we presented results about the variance, a 2 = 
(<5 2 ), and the skewness, 7 = (<5 3 )- 

It has been argued on several grounds (e.g. Coles & 
Frenk 1991) that the relation 



with 53 nearly independent of scale, should describe the 
clustering of cosmic structures. Although at small scales, be- 
low a few Mpc, eq. (23) is predicted by models of non-linear 
gravitational clustering (e.g. Borgani 1995 and references 
therein), at the larger scales, sampled by galaxy clusters, it 
is expected to hold due to mildly non-linear evolution as 
well as by the bias relating the cluster and DM distribu- 
tions. The resulting 53 values at different R sm are reported 
in column 6 of Table 3 for both the simulations and the 
Abell/ACO sample. We note that (a) the reduced skewness 
53 is always independent of the scale with a good accuracy, 
and (b) it takes the same value 53 ~ 1.9 for all the mod- 
els, within statistical fluctuations, and consistent with the 
observational data (PV95). 

Accordingly, we conclude that, for the cluster distribu- 
tion, only the amplitude of clustering, and not its nature, 
depends on the initial power-spectrum. This suggests that 
both the lognormal pdf shape and the 53 value observed for 
real data are natural consequences of high-peak biasing and 
possibly of the random-phase assumption of the primordial 
density field. 

5.3 The (3— parameter of clusters 

An interesting quantity, which relates the large-scale clus- 
tering to the linear peculiar velocity field of clusters is the 
/3-parameter, defined as 

A. = ^, (24) 

where /(f2 ) — 6 to a good accuracy (e.g. Peebles 1993) 
is the linear velocity factor and 



Q~cl 
0~DM 



(25) 



is the biasing parameter, here defined as the ratio between 
the r.m.s. fluctuations for the cluster and the DM density 
fields. 

For the simulations we directly estimate f3 c i, since we 
know a priori Q and we can compute (Jdm for each model 
according to 

-1 1/2 

dkk 2 P(k)W Rtm (k) 



0~DM{Rsm) 



(26) 



5 3 (a 2 ) 2 



(23) 



where Wji, m (fc) = exp( — k 2 R 2 m /2) is the Fourier transform 
of the window function of eq. (17). 

^From the observational size, f3 c i can be estimated by 
comparing the linear velocity induced by the cluster dis- 
tribution on the Local Group (LG) with the LG veloc- 
ity as measured from the CMB temperature dipole. Re- 
sults of these analyses (PV91; Scaramella, Vettolani & 
Zamorani 1991; Scaramella 1995; Branchini & Plionis 1995; 
Tini Brunozzi et al. 1995; Plionis 1995) consistently indicate 
that f3 c i = 0.20 ±0.05. From their power-spectrum analysis, 
Jing & Valdarnini (1993) and Peacock & Dodds (1994) ver- 
ified that at least the relative biasing between clusters and 
optically as well as infrared selected galaxies is independent 
of the scale to a quite good accuracy. In particular, Peacock 
& Dodds (1994) found that b c i/biRAS — 4.5, for the rela- 
tive biasing between clusters and IRAS galaxies. Therefore, 
if Piras = 1.0 ± 0.2 (Peacock & Dodds 1994), it turns out 
that f3 c i = 0.22 ± 0.05, in agreement with the other value. 
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In Figure 10 we plot the f3 c i parameters for the cluster 
simulations at different i? sm . For sake of clarity, we plot 
lcr error bars, due to cosmic variance, only for LOWH. The 
heavy-dashed horizontal lines delineates the lcr band for the 
observational result, f3 c i = 020 ±0.05. Two main conclusions 
can be drawn from this figure. 

(i) Independent of the model considered, f3 c i is fairly con- 
stant over the whole scale range, thus implying a nearly 
scale-independence also for b c i. This result is quite dif- 
ferent from the analogous one presented in Paper I, 
which showed a decreasing trend for b c i at small scales. 
This confirms how important is the increased resolu- 
tion and the optimization of the ZA for the reliability 
of our simulations. ^From one hand, such a linearity 
of the biasing is a rather remarkable result, since both 
the evolution of the density field and the selection of 
clusters as high-density peaks represent definitely non- 
linear transformations of the initial fluctuations. From 
the other hand, this result supports the usual assump- 
tion of b c i = const, used to infer the shape of the pri- 
mordial power-spectrum from that of clusters. 

(ii) Within the models with Q = 1, the TCDM is the only 
one that is only marginally consistent with the observa- 
tional f3 c i value, while all the other three models are per- 
fectly consistent. On the other hand, both the ACDM 
models give too small values of f3 c i. For ACDMi, the 
small f3 c i is due to the low spectrum normalization, 
which corresponds to a large b c i. Note, however, that 
for ACDM2 the f3 c i parameter is underestimated, due 
to the overestimate of the cluster correlations. In fact, 
if we would allow for ACDM2 clusters to have the same 
clustering as ACDMi clusters, then the f3 c i parameter 
would increase from f3 c i ~ 0.11 to f3 c i ~ 0.18, which 
is within the observational uncertainties. This indicates 
that, for a ACDM model, a normalization as high as 
as = 1.3 is required not only to match COBE data and 
cluster abundances, but also to generate a high enough 
f3 c i or, equivalently, to generate large enough peculiar 
velocities. 



6 CONCLUSIONS 

In this paper we have compared results of a statistical anal- 
yses of our simulated cluster distributions and of a com- 
bined sample of Abell/ACO clusters. Our cluster simula- 
tions, which are based on the truncated Zel'dovich approx- 
imation (TZA), are extremely cheap computationally; each 
realization takes about 6 minutes of CPU on a HP735/125 
workstation. This has allowed us to run a large number of 
realizations (50) for each model, so as to properly estimate 
the cosmic variance. 

In order to check the reliability of the TZA for clus- 
ter simulations, we compared them with analogous results 
based on PM N-body simulations, having the same initial 
conditions, mass resolution and cluster identification crite- 
ria as the TZA simulations. We find that, as long as cluster- 
ing is only mildly non-linear on the cluster mass scale (i.e. 
08 ~ 1), the TZA reproduces N-body results with remark- 
able accuracy. The agreement between the corresponding 
cluster distributions occurs not only in a statistical sense; 



Table 4. Summarizing scheme of the tests applied to the consid- 
ered DM models. The null score means that the test rejects the 
model at least at the 2a level. 

Model N(> M) J 3 pdf S 3 /3 c i 



SCDM 
TCDM 
LOWH 
CHDM 
ACDM 









1 1 



1 1 

1 1 

1 1 

1 1 
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cluster positions for TZA and PM simulations are shown to 
be extremely close point-by-point. Furthermore, comparing 
the resulting two-point correlation function for our CHDM 
clusters with that obtained by Klypin & Rhee (1994) from 
their higher resolution PM simulations of the same model, 
we find a extremely good agreement, even down to quite 
small scales (~ 7 ft. -1 Mpc), at which one can doubt the 
validity of the TZA. The reason for this remarkable suc- 
cess of the TZA in reproducing N-body results lies in its 
ability to account for non-local (i.e. long-wavelength) ef- 
fects when moving particles from their initial (linear) posi- 
tions to their correct evolved ones (cf. Pauls & Melott 1994; 
Sathyaprakash et al. 1994). 

However, if large spectrum normalizations are consid- 
ered (<J8^ 1), non-linear gravitational effects, like infall and 
merging of structures starts playing a role in N-body sim- 
ulations, while they are not accounted for in the TZA. As 
a result, the cluster clustering is stable for N-body simu- 
lations, while it increases with as in TZA simulations. We 
also verified that no substantial changes are found by going 
to the second-order in Lagrangian perturbative theory. This 
suggests the non-perturbative nature of effects like merging 
or infall. 

We applied applied five different tests to the considered 
DM models: (a) cluster abundances; (b) cluster correlations 
(Js(i?)); (c) shape of the cluster probability density function; 
(d) reduced skewness 53 ; (e) determination of the /3 C ; param- 
eter. In Table 6 we summarize the results. Failures, indicated 
with a '0', represent results incompatible with observations 
at the 2a level. For the ACDM model we report the results of 
the COBE-normalized version (ACDM2; as = 1.3) only for 
cluster abundances. Results concerning the clustering analy- 
ses are for the low-nomalization version (ACDMi ; as = 0.8), 
which are expected to be similar to those at as = 1.3 on the 
basis of stable clustering arguments, suggested by N-body 
simulations. 

As for the clustering analysis, we have analyzed both, 
the discrete cluster distribution using the J3 statistic, and 
the smoothed one using the pdf statistic. These analyses pro- 
vide stringent constraints on the DM models we have consid- 
ered and the only models that pass all the tests are the Cold 
+ Hot DM scenario and the low-density flat CDM model 
with the lower normalization [as = 0.8). Standard, Tilted 
and Low-Ho CDM versions do not account for the large- 
scale clustering of the real cluster sample. The low-density 
model with the larger normalization [as = 1.3) produces 
much stronger clustering than real clusters on all scales. 
However, in this case the TZA is expected to overestimate 
the cluster correlations. If, based on the clustering stability 
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suggested by N-body simulations, we accept that results at 
<78 = 0.8 are representative of those at as = 1.3, we could 
conclude that ACDM agrees quite well with the data for the 
clustering tests we have used. 

We find that the shape of the probability density func- 
tion (pdf) is not well approximated by a Gaussian distri- 
bution, even at the largest smoothing scale considered. The 
pdf for both simulations and Abell/ACO clusters is always 
much better reproduced by the lognormal model than by 
the Zel'dovich prediction, despite the fact that the ZA gov- 
erns the underlying dynamics. This shows that, at least for 
clusters, the lognormal shape of the pdf does not occur by 
chance, for a limited set of initial conditions, as argued to 
happen for the galaxy distribution (Bernardeau & Kofman 
1994). Instead it is much more likely to be related to the fact 
that clusters trace the high density peaks of the underlying 
matter field and, perhaps, to the initial random-phase as- 
sumption. Using a % 2- test to compare the shape of the pdf 
for both data and simulations, we find the best models to 
be CHDM and ACDM, consistent with the analysis of J3. 

Although the variance and skewness of the smoothed 
cluster pdf are powerful discriminators of different models 
(see Paper II), the reduced skewness, 53 = 7/cr 4 , turns out 
to be independent of the initial spectrum. We always find 
53 ~ 1.9, almost independent of the scale, and consistent 
with the observational results (Plionis & Valdarnini 1995). 
One could be tempted to conclude that such a value of 53 
is naturally produced by the high-peak selection of clusters, 
probably combined with the Gaussian nature of the initial 
fluctuations. Whether the analysis of the reduced skewness 
represents a test of Gaussian vs. non-Gaussian initial con- 
ditions remains to be seen (see also Coles et al. 1993). 

By computing the scale-dependence of the parameter 
fid = fJ?' 6 /&ci, we have verified that the linear biasing pre- 
scription used to relate cluster and DM distributions is al- 
ways satisfied to a good precision. The resulting value of f3 c i 
depends on the details of the model. Between the Q = 1 
models, the only one which can be marginally excluded is 
the TCDM, whose low normalization [as = 0.5) gives rise 
to a rather low f3 c i. Both the two low-density models have a 
quite low /3-parameter, f3 c i ~ 0.1. However, for ACDM2 this 
ought to be an underestimate, due to the overestimate of the 
clustering. On the other hand, taking the cluster correlations 
measured for ACDMi and rescaling to the normalization of 
ACDM2 , would give a result perfectly consistent with the 
observational constraint f3 c i = 0.20 ± 0.05. 

As a complementary indication of the success or not of 
the DM models considered, we compared the predicted clus- 
ter abundances with available observational results using the 
Press & Schechter (1974) approach. The only models that 
produce adequate abundances are CHDM, the low-density 
flat CDM model with <r 8 = 1.3 (ACDM 2 ) and low-H CDM. 
Given the uncertainties in both theory and observations, 
however, we take these results to be indicative rather than 
definitive. 

The overall picture emerging from these studies is that 
the large-scale cluster distribution places stringent con- 
straints on models of structure formation. Combining all 
the above results on clustering and abundances, it emerges 
that the only models to survive all the constraints we have 
considered here are CHDM and ACDM (with a marginal 
preference for the former). 



One can also ask whether reasonable modifications of 
the parameters in the DM models considered above (i.e. the 
Cold + Hot DM mixture, the values of Q and Qa values, 
the primordial spectral index n, the Hubble constant, etc.) 
could lead to significant changes in the resulting cluster dis- 
tribution. As an example, it has been shown that slightly 
decreasing the relative fraction of the hot component in the 
CHDM model bring it in a much better agreement with 
the high-redshift detection of collapsed structures (Klypin, 
Nolthenius & Primack 1995), while whether Qhot is shared 
between one or more neutrino species could also play a role. 
Verifying which of the CHDM variants produces an accept- 
able cluster distribution is not a difficult task, thanks to the 
low computational cost of our simulations. 

We believe that, although keeping in mind the limita- 
tions of the TZA, in future investigations of specific DM 
models, the optimal strategy would be, first of all, to run 
optimal TZA simulations in order to assess the model on 
large scales. Only after that one should decide whether a 
model is worth exploring at smaller scales by means of high- 
resolution, computationally expensive N-body simulations. 
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Figure Captions 

Figure 1. The average number of streams per Eulerian 
point, N s , as a function of the Gaussian filtering scale Rf 
for the six different models. The dotted horizontal line, 
N s = 1.1, delineates the quasi single-stream regime. 
Figure 2. Comparison between TZA and PM cluster simu- 
lations, at <78 = 0.67 (left panels) and as = 1 (right panels). 
In the upper panels the distribution of clusters obtained us- 
ing the two simulation techniques are compared. Open cir- 
cles are the TZA cluster positions (for Rf = 3 ft. -1 Mpc) and 
filled dots give the PM cluster positions. The lower panels 
compare the cluster count-in-cell variance for PM simula- 
tions and for TZA simulations, for three different values of 
the filtering radius Rf. 

Figure 3. Comparison between the two-point correlation 
function for our CHDM cluster simulations (open dots) and 
for the PM N-body results from Klypin & Rhee (1994) based 
on the same initial spectrum (filled dots). Their results are 
obtained with PM simulations having 256 3 grid points on 
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a box of side 200 ft. -1 Mpc. Our results are the average over 
50 realizations, and errors are la scatter over this ensemble. 
The KR94 results are an average over 2 realizations and the 
error bars are quasi-Poissonian estimates. 
Figure 4. The two-point cluster correlation functions for 
TZA simulations based on the first-order (filled circles) and 
second-order (open circles) Lagrangian theory, at as = 0.67 
(left panel) and as = 1 (right panel). The plotted £(r) values 
are the average over 10 realizations and the corresponding 
error bars are the la scatter over this ensemble. 
Figure 5. The cluster distribution (heavy dots) for the 
CHDM model superimposed on the DM particle distribution 
in a slice 10 h -1 Mpc thick for a box of side 640 h -1 Mpc. It 
is interesting to note that clusters are strongly correlated 
with the intersection of filaments of the DM distribution. 
Figure 6. The abundances of clusters with mass M > 4.2 x 
10 h~ Mq for the different DM models, as predicted by 
the Press & Schechter (1974) formalism, as a function of 
the critical density contrast S c . The horizontal lines are the 
observational results by White et al. (1993; dotted line) and 
by Biviano et al. (1993; dashed line) for clusters with mass 
larger than the above value. 

Figure 7. The Js(i?) integral as a function of the scale i? 
for the simulated (open circles) and real Abell/ACO (filled 
circles) cluster distributions. Error bars are plotted only for 
the simulations and correspond to la scatter over the en- 
semble of 50 realizations. 

Figure 8. Comparison between the pdfs for real (filled cir- 
cles) and simulated (open circles) cluster distributions. Re- 
sults only for the SCDM and CHDM are shown at R sm = 20 
and 40 ft. -1 Mpc. For the simulations, the analysis is realized 
in redshift space and the error bars correspond to cosmic 
r.m.s. scatter. 

Figure 9. Comparison between the pdfs of simulated 
cluster distributions and the theoretical models. We plot 
only results for SCDM and CHDM models at R sm = 
20 and 40 ft. -1 Mpc. Solid, long-dashed and short-dashed 
curves correspond to the lognormal, Zel'dovich and Gaussian 
model, respectively. Error bars are cosmic r.m.s. scatter. 
Figure 10. Scale dependence of the cluster /3-parameter 
for the six different models. For reasons of clarity, we plot 
the corresponding cosmic r.m.s. scatter only for the LOWH 
model. Similar uncertainties hold also for the other models. 
The heavy dashed horizontal lines show the range indicated 
by observational results (see text), f3 c i = 0.20 ± 0.05. 
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